{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definition of the Aerosol Optical Depth in terms of Angstrom Exponents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The authors of RRTM_SW use the following equation to determine the aerosol optical depth (AOD) at a given wavelength based on the optical depth at a reference wavelength (1 micron in RRTM_SW):\n",
    "\n",
    "$$ AOD(\\lambda) = \\frac{ AOD(\\lambda_{ref}) * (a_1 + a_2 * [\\frac{\\lambda}{\\lambda_{ref}}]) } {[(a_1 + a_2 -1) + [\\frac{\\lambda}{\\lambda_{ref}}]^{a_0} } $$\n",
    "\n",
    "Note that this formula contains three \"fit parameters\", $a_0$, $a_1$, and $a_2$, which allows one to describe a 2nd-order (parabolic) dependence of the AOD on wavelength.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "AERPAR      = np.array([2.184, 1.00, 0.00])\n",
    "wavelength  = np.arange(0.2, 3, 0.01)\n",
    "wavelength1 = 1.0\n",
    "#AOD1        = 0.046  # Total AOD = 0.25\n",
    "#AOD1        = 0.092  # Total AOD = 0.5\n",
    "AOD1        = 0.184  # Total AOD = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "specAOD = AOD1 * (AERPAR[1] + AERPAR[2] * (wavelength/wavelength1)) / ((AERPAR[1] + AERPAR[2] - 1) + (wavelength/wavelength1)**AERPAR[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1211af9b0>"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtkAAAHgCAYAAABw0HFmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3hVVcLF4bXTSSEhhCSkEAJBIAQCJPQOSlGKioqC2EUdu+M444wz6ozz2XsHuyLYlagD0qJIk957DyWhBkJNyP7+IDIMQgyQm3PL732e88g9995kwTHJYrPP3sZaKwAAAACVx8/pAAAAAIC3oWQDAAAAlYySDQAAAFQySjYAAABQySjZAAAAQCWjZAMAAACVLMDpAK4QExNj69at63QMn7V//36FhYU5HQNngGvmebhmnodr5nm4Zp6nqq/ZnDlzdlhra53qOa8s2XXr1tXs2bOdjuGzcnNz1bVrV6dj4AxwzTwP18zzcM08D9fM81T1NTPGbDjdc141XcQY088YM7ywsNDpKAAAAPBhXlWyrbU51tphkZGRTkcBAACAD/Oqkg0AAAC4A6+ckw0AAICKKy4uVl5eng4dOuR0lHMSGRmpZcuWVfrHDQkJUVJSkgIDAyv8Hko2AACAj8vLy1NERITq1q0rY4zTcc7avn37FBERUakf01qrnTt3Ki8vT6mpqRV+H9NFAAAAfNyhQ4dUs2ZNjy7YrmKMUc2aNc94lJ+SDQAAAAp2Oc7mz4aSDQAAAMeFh4c79rl79+6tqKgo9e3bt9I+JiUbAAAAPqOkpOQ35/70pz/pww8/rNTP41Ulm81oAAAAvEdOTo7atGmjFi1a6Pzzz1d+fr5KS0vVoEEDbd++XZJUWlqqtLQ07dixQzt27NDAgQPVqlUrtWrVSlOnTpUkPfLIIxo2bJh69uypa6655jefp0ePHpV+w6RXrS5irc2RlJOdnX2z01kAAAA80aM5S7R0y95K/ZjpCdX1cL8mZ/y+jh07asaMGTLG6K233tJTTz2lZ599VldffbVGjhype+65RxMmTFBmZqZiYmJ022236d5771XHjh21ceNG9erV6/iSfnPmzNHPP/+satWqVerv7XS8qmQDAADAe+Tl5WnQoEHaunWrjhw5cnwJvRtuuEEDBgzQPffco3feeUfXX3+9JCk3N1erVq06/v69e/dq3759kqT+/ftXWcGWKNkAAAA4wdmMOLvKnXfeqfvuu0/9+/dXbm6uHnnkEUlScnKy4uLiNGnSJM2cOVMjR46UdGzqyPTp009ZpsPCwqoyunfNyQYAAID3KCwsVGJioiTp/fff/5/nbrrpJl199dW64oor5O/vL0nq3r27XnnlleOvmT9/ftWFPQklGwAAAI47cOCAkpKSjh/PPfecHnnkEV1++eXq1KmTYmJi/uf1/fv3V1FR0fGpIpL09NNPa/bs2WrWrJnS09P1xhtvVOhzd+rUSZdffrkmTpyopKQkjRs37px/P0wXAQAAgONKS0tPeX7AgAGnPL9gwQJlZmaqUaNGx8/VrFlTn3zyyW9e++s0k9OZMmVKxYNWECPZlaS01Gpr4UGnYwAAAHi9J554QgMHDtTjjz/udJTTomRXkrs/ma8rh8+QtdbpKAAAAF7tL3/5izZs2KCOHTs6HeW0KNmVpE1qtDbsPKBVBUVORwEAAIDDKNmV5IL0OEnSD0u2OZwEAADgzPGv8ad3Nn82XlWyndxWPa56iJonR+mHpflV/rkBAADORUhIiHbu3EnRPgVrrXbu3KmQkJAzep9XrS7i9LbqvZrE68mxy7Vlz0ElRFXdjkIAAADnIikpSXl5edq+fbvTUc7JoUOHzrgMV0RISIiSkpLO6D1eVbKd1rNJnJ4cu1zjl+br2vZ1nY4DAABQIYGBgce3LPdkubm5atGihdMxJHnZdBGn1a8Vrvq1wvTDUuZlAwAA+DJKdiXr2SReM9buUuGBYqejAAAAwCGU7ErWMz1OR0utJq3gBkgAAABfRcmuZJlJUYqNCNYPSyjZAAAAvoqSXcn8/Ix6NonTjyu361DxUafjAAAAwAGUbBfomR6vA0eOaurqHU5HAQAAgAMo2S7Qtl5NRQQHMGUEAADAR1GyXSAowE/dGsVqwrJ8HS1l5yQAAABfQ8l2kZ5N4rRz/xHNWr/L6SgAAACoYl5Vso0x/YwxwwsLC52Oom4NYxUS6KfvFm51OgoAAACqmFeVbGttjrV2WGRkpNNRFBYcoB6N4/T9oq0qOVrqdBwAAABUIa8q2e6mX7ME7dx/RNPX7nQ6CgAAAKoQJduFujaspfDgAOUs2OJ0FAAAAFQhSrYLhQT6q2eTOI1dvE2HS9iYBgAAwFdQsl2sX2aC9h4q0ZSVbEwDAADgKyjZLtYxLUZRoYHKWciUEQAAAF9ByXaxQH8/9cmI1/il+Tp4hCkjAAAAvoCSXQX6NUvQgSNHNXlFgdNRAAAAUAUo2VWgTb2aigkPZpURAAAAH0HJrgL+fkZ9m9XWpOUF2neo2Ok4AAAAcDFKdhXpl1lbh0tKNWFZvtNRAAAA4GKU7CrSIrmGEqOq6at5TBkBAADwdpTsKuLnZ3RJi0T9vGq78vcecjoOAAAAXIiSXYUGZiWp1EpfzdvsdBQAAAC4kFeVbGNMP2PM8MLCQqejnFJqTJiyUmro8zl5stY6HQcAAAAu4lUl21qbY60dFhkZ6XSU07osK0mrC4q0MM89/yIAAACAc+dVJdsTXNSstoID/PT5nDynowAAAMBFKNlVrHpIoHo1ideYBVt0uIRt1gEAALwRJdsBA7OSVHiwWBOXsc06AACAN6JkO6BjWoziqgfrC6aMAAAAeCVKtgP8/YwuaZGk3JXbtX3fYafjAAAAoJJRsh1yWVaijpZafTOfNbMBAAC8DSXbIWmxEcpMjmLNbAAAAC9EyXbQ5VlJWr5tn+Zv2uN0FAAAAFQiSraDBjRPUGiQvz6eudHpKAAAAKhElGwHRYQEakDzBOUs3KLCg8VOxwEAAEAloWQ7bHDrFB0qLtXX87gBEgAAwFtQsh3WNClSTRMj9fHMjdwACQAA4CUo2W5gcJs6WpG/T3M37nY6CgAAACoBJdsN9M9MUHhwgEZyAyQAAIBXoGS7gbDgAF3cIkHfLtyqPQeOOB0HAAAA54iS7SYGt07RkZJSfTGXGyABAAA8HSXbTaQnVFfz5Ch9PHMDN0ACAAB4OK8q2caYfsaY4YWFhU5HOSuD29TRmu37NXPdLqejAAAA4Bx4Vcm21uZYa4dFRkY6HeWs9GuWoKjQQL0/bb3TUQAAAHAOvKpke7pqQf66slUdjVuyTXm7DzgdBwAAAGeJku1mhrZLkTFGH07f4HQUAAAAnCVKtptJjKqmXk3iNOqXjTpwpMTpOAAAADgLlGw3dH2HVO09VKKv5rGcHwAAgCeiZLuh7JQaykisrvemrmc5PwAAAA9EyXZDxhhd1z5VqwqK9PPqHU7HAQAAwBmiZLupfpm1FRMepHenrnc6CgAAAM4QJdtNBQf4a3CbFE1aXqB1O/Y7HQcAAABngJLtxq5uW0eB/obNaQAAADwMJduNxUaEqF9mgj6ZtUm79x9xOg4AAAAqiJLt5oZ1rqeDxUf14Qw2pwEAAPAUlGw31yi+uro1rKX3pq3XoeKjTscBAABABVCyPcCtXepr1/4j+mz2JqejAAAAoAIo2R6gdWq0midHacSUdSo5Wup0HAAAAPwOSrYHMMbo1i71tXHXAf1n8Tan4wAAAOB3ULI9xAXpcaoXE6Y3f1rDVusAAABujpLtIfz9jG7uXE+LN+/V1NU7nY4DAACAclCyPcglLRJVKyJYb/60xukoAAAAKAcl24OEBPrr+g51NWXVDi3YtMfpOAAAADgNSraHGdo2RZHVAvXypFVORwEAAMBpULI9TERIoG7smKoJywq0eHOh03EAAABwCpRsD3Rt+7qKCAlgNBsAAMBNUbI9UGS1QF3fIVXjluRr2da9TscBAADASSjZHuqGDnUVHhygVyatdjoKAAAATkLJ9lBRoUG6tn2Kvl+8VSvz9zkdBwAAACfwqpJtjOlnjBleWOgbNwTe2LGeqgX6M5oNAADgZryqZFtrc6y1wyIjI52OUiWiw4I0tF2KchZu0eqCIqfjAAAAoIxXlWxfdHOnegoJ8NeLE1lpBAAAwF1Qsj1cTHiwru9QVzkLtmjpFlYaAQAAcAeUbC9wS+f6iggJ0HPjVzgdBQAAAKJke4XI0EDd2qW+Jiwr0NyNu52OAwAA4PMo2V7iuvZ1FRMepGfGMZoNAADgNEq2lwgLDtAfuqZp2pqdmrp6h9NxAAAAfBol24sMblNHCZEhemrcCllrnY4DAADgsyjZXiQk0F93n99ACzbt0YRlBU7HAQAA8FmUbC8zsGWSUmPC9PS45So5Wup0HAAAAJ9EyfYyAf5++lOvhlqZX6TP5+Q5HQcAAMAnUbK9UJ+MeLWsE6Xnxq/UgSMlTscBAADwOZRsL2SM0d8uaqyCfYc14qd1TscBAADwOZRsL5WVEq0+GfF686c1Kth3yOk4AAAAPoWS7cX+3LuRjpSU6vnxq5yOAgAA4FMo2V6sbkyYrm6bok9mbdSq/H1OxwEAAPAZlGwvd1ePBgoLDtDj/1nudBQAAACfQcn2ctFhQbq9W5omLS/QlFXbnY4DAADgEyjZPuC69nVVJzpU/8xZqmI2qAEAAHA5SrYPCAn010MXNdaqgiJ9NGOD03EAAAC8HiXbR1yQHqdODWL0/PiV2ll02Ok4AAAAXo2S7SOMMfpH33TtP3JUz45f6XQcAAAAr0bJ9iEN4iJ0TbsUjfplo5ZsKXQ6DgAAgNeiZPuYe84/TzVCg/TomKWy1jodBwAAwCtRsn1MZLVA3d+zoX5Zv0s5C7c6HQcAAMArUbJ90KBWycpIrK7Hvl2qfYeKnY4DAADgdSjZPsjfz+ixi5tqe9FhPT9+ldNxAAAAvA4l20c1T47S4NZ19N60ddwECQAAUMko2T7sgV6NVCM0SH//erFKS7kJEgAAoLJQsn1YZGig/nphY83duEefzt7kdBwAAACvQcn2cZe2TFTr1Gg9MXa5du0/4nQcAAAAr0DJ9nHGGD12cYaKDpXo8e+XOR0HAADAK1CyofPiInRz53r6bE6epq3e4XQcAAAAj0fJhiTp7h4NlFIzVA9+tUiHio86HQcAAMCjUbIhSQoJ9NfjlzbVhp0H9MIE1s4GAAA4F5RsHNe+fowGZSdrxJS1WryZtbMBAADOFiUb/+OvFzZWdFiQ/vzFQpUcLXU6DgAAgEeiZON/RIYG6p/9m2jJlr16++d1TscBAADwSJRs/EbvjHj1TI/Tc+NXat2O/U7HAQAA8DiUbPyGMUb/ujhDwQF++tNnC3SULdcBAADOCCUbpxRXPUSP9G+i2Rt2692pTBsBAAA4E5RsnNYlLRJ1fuM4PT1uhdZsL3I6DgAAgMdw+5JtjAkzxrxvjBlhjBnidB5fYozR/12aoWpB/rqfaSMAAAAV5kjJNsa8Y4wpMMYsPul8b2PMCmPMamPMX8pOXyrpc2vtzZL6V3lYHxcbEaJH+zfRvI179NaUtU7HAQAA8AhOjWS/J6n3iSeMMf6SXpXUR1K6pKuMMemSkiRtKnsZ+307oH9mgno3idez41dqVf4+p+MAAAC4PUdKtrX2J0m7TjrdWtJqa+1aa+0RSaMlDZCUp2NFW/KA6S3eyBijxy7JUHhwgO79dL6OlLBJDQAAQHmMtc7MszXG1JX0rbU2o+zxZZJ6W2tvKns8VFIbSX+W9IqkQ5J+ttaOPM3HGyZpmCTFxcVljR492tW/BZ8zJ79EL887rL71AnXZeUGnfV1RUZHCw8OrMBnOFdfM83DNPA/XzPNwzTxPVV+zbt26zbHWZp/quYAqS/H7zCnOWWvtfknX/96brbXDJQ2XpOzsbNu1a9fKTQd1lZTvv0Cfz8nTtT2z1apu9Clfl5ubK/78PQvXzPNwzTwP18zzcM08jztdM3eafpEnKfmEx0mStjiUBafxj35NlFQjVPd+Ml/7DhU7HQcAAMAtuVPJniWpgTEm1RgTJOlKSWMczoSThAcH6PlBzbVlz0E9Mmap03EAAADcklNL+I2SNF1SQ2NMnjHmRmttiaQ7JI2TtEzSp9baJU7kQ/myUmro9m5p+mJunr5ftNXpOAAAAG7HkTnZ1tqrTnP+e0nfV3EcnIW7ejTQjyu368EvF6l5cpQSoqo5HQkAAMBtuNN0EXiQQH8/vXhlC5UcLdVdo+ap5CjL+gEAAPzKq0q2MaafMWZ4YWGh01F8QmpMmP59SVPN3rBbL0xY5XQcAAAAt+FVJdtam2OtHRYZGel0FJ9xcYtEXZ6VpFdzV+vnVTucjgMAAOAWvKpkwxmPDmii+rXCdc8n87V932Gn4wAAADiOko1zFhoUoFcGt9C+Q8W679P5KnVoF1EAAAB3QclGpWgUX10P92uiKat26Pt1bFIDAAB8GyUbleaq1sm6qFltfbmqWLPX73I6DgAAgGMo2ag0xhg9fmlT1QwxumvUPO05cMTpSAAAAI6gZKNSVQ8J1G3Ng7W96LDu+3SBSkuZnw0AAHyPV5Vs1sl2D/Ui/fX3vumatLxAr0xe7XQcAACAKudVJZt1st3H0LYpuqRFop6fsFKTVxQ4HQcAAKBKeVXJhvswxuj/LmmqhnERumf0fG3adcDpSAAAAFWGkg2XqRbkrzeHZslaq1s+nKNDxUedjgQAAFAlKNlwqZSaYXrhyuZaunWv/vbVYlk2qgEAAD6gQiXbGONvjEkwxtT59XB1MHiP7o3idHePBvpibp5GztzodBwAAACXC/i9Fxhj7pT0sKR8SaVlp62kZi7MBS9zd48GWpC3R4/mLFF6QnW1rFPD6UgAAAAuU5GR7LslNbTWNrHWNi07KNg4I35+Ri8Maq74yBD94aO5Kth3yOlIAAAALlORkr1JEgtP45xFhQbpjauzVHiwmBshAQCAVzttyTbG3GeMuU/SWkm5xpgHfz1Xdt7tsBmN+2uSEKnnB2Vq3sY9+uuXi7gREgAAeKXyRrIjyo6NksZLCjrhXLjro505NqPxDL0zauuPF5ynL+dt1ps/rXU6DgAAQKU77Y2P1tpHJckYc7m19rMTnzPGXO7qYPBud3RP08qCIj05drnSaoXr/PQ4pyMBAABUmorMyX6wgueACjPG6OnLmqlpYqTuHj1PK7btczoSAABApSlvTnYfY8zLkhKNMS+dcLwnqaTKEsJrhQT6a/jQbIUFB+imD2Zp1/4jTkcCAACoFOWNZG+RNFvSIUlzTjjGSOrl+mjwBfGRIRp+TbYK9h7WbR/N0ZGS0t9/EwAAgJs7bcm21i6w1r4vKU3SKEnzJM2V9K21dncV5YMPaJ4cpacua6aZ63bpoa9ZcQQAAHi+393xUdIFkt6UtEaSkZRqjLnFWvsflyaDTxnQPFFrtu/XSxNXKblGqO7s0cDpSAAAAGetIiX7OUndrLWrJckYU1/Sd5Io2ahU957fQHm7DujZ8SuVFF1Nl7RIcjoSAADAWalIyS74tWCXWSupwEV54MOMMXpiYDNt23tID3y+UHERIWqfFuN0LAAAgDNWkSX8lhhjvjfGXGeMuVZSjqRZxphLjTGXujjfGWHHR88XFOCn16/OUmpMmG75aI5W5rO0HwAA8DwVKdkhkvIldZHUVdJ2SdGS+knq67JkZ4EdH71DZLVAvXt9a1UL9Nd17/yi/L2HnI4EAABwRn53uoi19vqqCAKcKDGqmt65rpWueHO6bnhvlj69pZ3CgisyuwkAAMB5vzuSbYw5zxgz0RizuOxxM2PMQ66PBl+XkRipV4e01PJt+3T7x3NVfJQ1tAEAgGeoyHSRETq2jXqxJFlrF0q60pWhgF91axirxy7OUO6K7frzFwtVWsoa2gAAwP1V5N/fQ621vxhjTjzHtuqoMle1rqPt+w7rufErFR0apL9d1Fgn/f8IAADgVipSsneUrY1tJckYc5mkrS5NBZzkzu5p2ll0WG/9vE41w4N1W9f6TkcCAAA4rYqU7NslDZfUyBizWdI6SUNcmgo4iTFGD/drot0HivXk2OWKDgvUoFZ1nI4FAABwShVZXWStpPONMWGS/Ky1LFwMR/j5GT1zeab2HCzWg18uUmS1IPXOiHc6FgAAwG+Ue+OjMaahMeZZY8x3kj6V9LAx5ryqiQb8VlCAn964uqUyk6N01+h5mr5mp9ORAAAAfuO0JdsY005SrqR9OjZdZISk/ZJyjTFtqyQdcAqhQQF697pWSokO1c0fzNbizezwCQAA3Et5I9n/kHSVtfYRa+031tqvrbUPS7pK0sNVEw84tajQIH14YxtFVgvUNe/8wvbrAADArZRXsutba3NPPmmt/VFSPZclOgfGmH7GmOGFhYxs+oL4yBCNvKmNAvyMhrw1U+t27Hc6EgAAgKTyS3Z5Q4Nu2WastTnW2mGRkZFOR0EVqRsTppE3tdHRUqshI2Zo064DTkcCAAAod3WRZGPMS6c4byQluigPcMYaxEXooxvb6Mrh0zX4rRn67Jb2io8McToWAADwYeWNZP9J0pxTHLMlPeD6aEDFpSdU14c3ttHu/cUa/NYMbd932OlIAADAh512JNta+35VBgHOVWZylN69vpWuefsXXf3WTI0e1lY1woKcjgUAAHxQuetkA56mVd1ovXVtttbt3K+h78xU4cFipyMBAAAfRMmG1+mQFqM3r87Sim37NPTtmSo8QNEGAABVi5INr9StUazeuDpLy7fu05C3Z2jPgSNORwIAAD7ktHOyjTEvS7Kne95ae5dLEgGVpEfjOL05NEu3fDhHQ96aqZE3tVFUKHO0AQCA65U3kj1bp15d5NcDcHvdGsVq+DVZWlVQpMEjZmr3fka0AQCA67G6CLxe14axGnFNtm7+YLYGl41oR7PqCAAAcKHfnZNtjKlljHnGGPO9MWbSr0dVhAMqS5fzaumta7K1dnuRBo+YoV2MaAMAABeqyI2PIyUtk5Qq6VFJ6yXNcmEmwCU6n1fr2PJ+O/brquEzVLDvkNORAACAl6pIya5prX1bUrG19kdr7Q2S2ro4F+ASnRrU0rvXtdKm3Qd0xRvTtXnPQacjAQAAL1SRkv3rIsNbjTEXGWNaSEpyYSbApdqnxejDG9to5/4juuKN6Vq3Y7/TkQAAgJepSMl+zBgTKemPku6X9Jake12a6iwZY/oZY4YXFhY6HQVuLiulhkbd3FYHi4/qijena8W2fU5HAgAAXuR3S7a19ltrbaG1drG1tpu1NstaO6Yqwp0pa22OtXZYZGSk01HgATISI/XJsLbyM9Kg4dO1MG+P05EAAICXqMjqIu8bY6JOeFzDGPOOa2MBVaNBXIQ+u6W9woMDNHjETP2ybpfTkQAAgBeoyHSRZtba40N81trdklq4LhJQterUDNVnt7ZTbPVgXfPOTOWuKHA6EgAA8HAVKdl+xpgavz4wxkSrnE1sAE9UO7KaPr2lnerFhOum92fr63mbnY4EAAA8WEVK9rOSphlj/mWM+ZekaZKecm0soOrFhAdr9C1tlV23hu75ZL7e/nmd05EAAICHqsiNjx9IGigpX1KBpEuttR+6OhjghOohgXrv+tbq3SRe//p2qZ4cu1zWWqdjAQAAD3Pakm2MqV7232hJ2yR9rGO7P24rOwd4pZBAf706pKUGt6mj13PX6M9fLFTJ0VKnYwEAAA9S3tzqjyX1lTRH0olDeabscT0X5gIc5e9n9O+LMxQTHqyXJq7Srv3FemVwC4UE+jsdDQAAeIDTlmxrbd+y/6ZWXRzAfRhjdN8F56lmWJAeyVmioW/P1IhrshUVGuR0NAAA4OYqsk72xIqcA7zVte3r6uWrWmjBpkJd+vo0bdp1wOlIAADAzZU3JzukbO51TNkGNNFlR11JCVUVEHAHfZsl6KOb2mhn0RFd8to0docEAADlKm8k+xYdm4/dSNLcsl/PkfSNpFddHw1wL61To/XFbe0UEuinQW/O0MRl+U5HAgAAbuq0Jdta+2LZfOz7rbWpJxyZ1tpXqjAj4DbSYiP05R/aKy02XDd/MFsjZ25wOhIAAHBDFdmM5i1jzH3GmC+NMV8YY+4xxoS4PBngpmIjQjR6WFt1bRirv321WE+OXa7SUtbSBgAA/1WRkv2+pCaSXpb0iqR0SWxGA58WFhyg4UOzNKRsLe07R83TwSNHnY4FAADcRHnrZP+qobU284THk40xC1wVCPAUAf5+euziDKXUDNXj/1muTbsPaMQ12Yqrzj/0AADg6yoykj3PGNP21wfGmDaSprouEuA5jDEa1rm+hg/N1uqCIg14ZaoWby50OhYAAHBYRUp2G0nTjDHrjTHrJU2X1MUYs8gYs9Cl6QAPcUF6nL64rb38/Ywue2Oa/rNoq9ORAACAgypSsntLSpXUpexIlXShjm253s910c6cMaafMWZ4YSEjiah6jWtX19e3d1Dj2tV128i5emXSKlnLDZEAAPiiipTsejpWpvtKqmet3XDi4dp4Z8Zam2OtHRYZGel0FPioWhHBGnVzW13cPEHP/LBS934yX4eKuSESAABfc9obH40xiZK+lHRIxzahMZKuMMY8KekSa+3mqokIeJaQQH89P6i50mLD9cwPK7Vx1wG9OTRbtSKCnY4GAACqSHkj2a9Iet1a28Vae5+19l5rbZey869VTTzAMxljdEf3Bnp9SEst3bpXF7/KDZEAAPiS8kp2urX2vZNPWms/0LGt1gH8jj5Na+uzW9qr1FoNfH2avpqX53QkAABQBcor2f6nOmmM8TvdcwB+q2lSpHLu7KjM5Cjd+8kC/TNnqYqPljodCwAAuFB5JTvHGDPCGBP264myX78h6XuXJwO8SEx4sEbe1EbXta+rd6au09C3Z2pn0WGnYwEAABcpr2Q/IKlQ0gZjzBxjzGxJ6yXtlXR/FWQDvEqgv58e6d9Ez16eqXkb96jfyz9rUR7ztAEA8EanLdnW2mJr7f2SkiVdJ+kGSSnW2vuttUeqKB/gdQZmJenzW9vLGKOBb0zTF3OYpw0AgLf53XWyrbUHrbWLrLULrbUHqiIU4O2aJjW32A0AACAASURBVEVqzB0dlFWnhv742QI9MmYJ87QBAPAiFdmMBoAL1AwP1oc3ttaNHVP13rT1GjJipvL3HnI6FgAAqASUbMBBAf5++nvfdL14ZXMt2lyoi16aommrdzgdCwAAnKPydnxsWd4brbVzKz8O4JsGNE9Ueu3qum3kXF399kzde/55ur1bmvz8jNPRAADAWThtyZb0bDnPWUndKzkL4NMaxEXom9s76K9fLdKz41dq1obdemFQc0WHBTkdDQAAnKHTlmxrbbeqDAJACgsO0AuDmqt1arQeHbNUF700Ra8MbqmslBpORwMAAGegQnOyjTEZxpgrjDHX/Hq4Ohjgq4wxGtImRV/+ob0C/I0GvTldb01ZK2ut09EAAEAF/W7JNsY8LOnlsqObpKck9XdxLsDnZSRG6ts7O6l7o1g99t0y3frRHBUeLHY6FgAAqICKjGRfJqmHpG3W2uslZUoKdmkqAJKkyGqBenNolh66qLEmLitQ35enaN7G3U7HAgAAv6MiJfugtbZUUokxprqkAkn1XBsLwK+MMbqpUz19cks7lZZKl78xXa/lrlZpKdNHAABwVxUp2bONMVGSRkiaI2mupF9cmgrAb2Sl1ND3d3dSrybxemrsCg19Z6YK2LwGAAC3VG7JNsYYSY9ba/dYa9+QdIGka8umjQCoYpHVAvXK4BZ64tKmmrNht3q/OEWTlxc4HQsAAJyk3JJtjy1n8PUJj9dbaxe6PBWA0zLG6MrWdfTtnR0VGxGs69+bpX99u1SHS446HQ0AAJSpyHSRGcaYVi5PAuCMpMVG6OvbO+jadil6++d1Gvj6NK3dXuR0LAAAoIqV7G6Sphtj1hhjFhpjFhljGM0G3EBIoL8eHZCh4UOzlLf7oPq+/LM+mbWRNbUBAHBYeduq/6qPy1MAOCc9m8SraVKk7v1kvv78xSJNXFagxy9tqprhrLYJAIATKjKS/Zi1dsOJh6THXB0MwJmpHVlNH9/UVn+7sLFyV2xXrxemaNLyfKdjAQDgkypSspuc+MAY4y8pyzVxAJwLPz+jmzvX05g7OygmPEg3vDdbf/tqkQ4cKXE6GgAAPuW0JdsY86AxZp+kZsaYvWXHPh3bjOabKkt4Bowx/YwxwwsLC52OAjiqUXx1fXNHBw3rXE8f/7JRF730s+Zv2uN0LAAAfMZpS7a19nFrbYSkp6211cuOCGttTWvtg1WYscKstTnW2mGRkZFORwEcFxzgr79e2Fgf39RWh4uPauDr0/TChJUqOVrqdDQAALxeRaaL/GKMOd5ajTFRxpiLXZgJQCVqV7+m/nNPZ/XPTNALE1Zp4BvTtbqApf4AAHClipTsh621x+dfWGv3SHrYdZEAVLbIaoF6flBzvTK4hTbs3K8LX5qiET+t1dFSlvoDAMAVKlKyT/Waiiz9B8DN9G2WoB/u7awu59XSv79fpivenM4GNgAAuEBFSvZsY8xzxpj6xph6xpjnJc1xdTAArhEbEaLhQ7P0wqDmWl1QpD4vTtG49cWMagMAUIkqUrLvlHRE0ieSPpV0UNLtrgwFwLWMMbq4RaLG39tZnRrEaNTyI7py+HSt37Hf6WgAAHiF3y3Z1tr91tq/SOpqrc221v7VWstPYsALxFYP0YhrsnVz0yCt2LZPvV/8Se9OXadSRrUBADgnv1uyjTHtjTFLJS0te5xpjHnN5ckAVAljjDokBuqHe7uoff0YPZqzVFcOn6E1zNUGAOCsVWS6yPOSeknaKUnW2gWSOrsyFICqFx8ZorevzdbTlzXT8m171efFKXp18moVs642AABnrCIlW9baTSedOuqCLAAcZozR5dnJmvDHLjq/cayeHrdC/V+ZqoV57BYJAMCZqEjJ3mSMaS/JGmOCjDH3S1rm4lwAHBQbEaLXhmTpzaFZ2ll0WBe/OlX/9/0yHTzC368BAKiIipTsW3VsNZFESZslNReriwA+oVeTeI2/r4sGtUrW8J/WqveLP2na6h1OxwIAwO1VZHWRHdbaIdbaOGttLWvt1dbanVURDoDzIqsF6vFLm2nUzW1lJA1+a6b+/PlCFR4odjoaAABuqyKri9QzxuQYY7YbYwqMMd8YY+pVRTgA7qNd/Zoae09n3dKlnj6fm6cez+Xq63mbZS3L/QEAcLKKTBf5WMc2oaktKUHSZ5JGuTIUAPcUEuivB/s01je3d1BijVDd88l8DXlrJsv9AQBwkoqUbGOt/dBaW1J2fCSJoSvAh2UkRurL29rrsYsztGhzofq8MEXPjV+pQ8XcGAkAgFSxkj3ZGPMXY0xdY0yKMeYBSd8ZY6KNMdGuDgjAPfn7GV3dNkUT/9hFfZrG66WJq9T7hZ80ZdV2p6MBAOC4ipTsQZJukTRZUq6k2yTdIGmOpNkuSwbAI8RGhOjFK1vooxvbyBijoW//ortGzVPBvkNORwMAwDEBv/cCa21qVQQB4Nk6NojRf+7upDd+XKPXJq/R5OUFeqB3Qw1ukyJ/P+N0PAAAqtRpR7KNMa2MMfEnPL6mbGWRl5gmAuBUQgL9dc/552ncvZ2VmRylv3+zRJe+xo6RAADfU950kTclHZEkY0xnSU9I+kBSoaThro8GwFOlxoTpwxtb68Urm2vznkMa8OpUPfjlQu0sOux0NAAAqkR5JdvfWrur7NeDJA231n5hrf27pDTXRwPgyYwxGtA8UZPu76IbO6Tqs9l56vZMrt6ftl4lR0udjgcAgEuVW7KNMb/O2e4hadIJz/3uXG4AkKTqIYF6qG+6xt7TSc2SovTwmCXq+/LPmrmWjWMBAN6rvJI9StKPxphvJB2UNEWSjDFpOjZlBAAqLC02Qh/e2FqvD2mpfYdKNGj4DN01ap62FbIKCQDA+5x2RNpa+29jzEQd2+nxB/vfvZP9JN1ZFeEAeBdjjPo0ra2uDWP1+o9r9MaPazRhWb7u6J6mGzumKjjA3+mIAABUinLXybbWzrDWfmWt3X/CuZXW2rmujwbAW1UL8td9F5ynifd1Uce0GD01doV6vzBFk5cXOB0NAIBKUZHNaADAJZKjQzX8mmx9cENrGSNd/94sXffuL1pdsM/paAAAnBNKNgDHdT6vlsbe3VkPXdRYczbsVq8Xpugf3yzWrv1HnI4GAMBZoWQDcAtBAX66qVM9/finbhrSpo5GztyoLk9P1oif1upwyVGn4wEAcEYo2QDcSnRYkP45IENj7+6krJQa+vf3y9Tz+Z80dvFW/ff+awAA3BslG4BbahAXofeub633b2it4AA/3frRXF05fIYWb2YFUQCA+6NkA3BrXc6rpe/v6qTHLs7Q6oIi9XvlZ/3x0wXK38v62gAA90XJBuD2Avz9dHXbFE3+U1cN61xPOQu2qOvTuXruhxUqOlzidDwAAH6Dkg3AY1QPCdSDfRprwn1d1L1xrF6atFpdnpqs96et15GSUqfjAQBwHCUbgMepUzNUrw5uqW9u76AGceF6eMwSXfD8j8pZsIWbIwEAboGSDcBjZSZHadTNbfXuda0UEuCvO0fN04BXp2ramh1ORwMA+DhKNgCPZoxRt0ax+v7uTnrm8kzt2HdYg0fM1HXv/qLl2/Y6HQ8A4KMo2QC8gr+f0WVZSZp0f1c92KeR5m7YrT4vTtEfP12gzXsOOh0PAOBjKNkAvEpIoL9u6VJfPz3QTTd3qqechVvU7Zlc/evbpdpZdNjpeAAAH0HJBuCVokKD9NcLG2vSH7uof2aC3p26Tp2emqxnf1ihwoPFTscDAHg5SjYAr5ZUI1TPXJ6pH+7tom6NYvXypNXq/NRkvZa7WgeOsMY2AMA1KNkAfEJabLheHdxS397ZUVkpNfTU2BXq/FSu3p26TodLjjodDwDgZdy+ZBtj6hlj3jbGfO50FgCeLyMxUu9c10pf3NZOabFhejRnqbo/86M+mbVRJUfZ0AYAUDlcWrKNMe8YYwqMMYtPOt/bGLPCGLPaGPOX8j6GtXattfZGV+YE4HuyUqI16ua2+ujGNooJD9Kfv1ikns//pDELtqi0lA1tAADnxtUj2e9J6n3iCWOMv6RXJfWRlC7pKmNMujGmqTHm25OOWBfnA+DDjDHq2CBGX9/eQcOHZinQ3093jZqnXi/8pG8XUrYBAGcvwJUf3Fr7kzGm7kmnW0taba1dK0nGmNGSBlhrH5fU15V5AOBUjDHq2SRePRrH6ftFW/XSxFW64+N5ahC7Snef30AXZtSWn59xOiYAwIMYa107UlNWsr+11maUPb5MUm9r7U1lj4dKamOtveM0768p6d+SLpD0VlkZP9XrhkkaJklxcXFZo0ePruTfCSqqqKhI4eHhTsfAGeCa/a9SazVr21F9s+aIthRZJYQbXVw/SNnx/vIz7lG2uWaeh2vmebhmnqeqr1m3bt3mWGuzT/WcS0eyT+NUP6FO2/SttTsl3fp7H9RaO1zScEnKzs62Xbt2Pdt8OEe5ubniz9+zcM1+q7uk+0vt8ZHt1xYUqcHWcLcZ2eaaeR6umefhmnked7pmTqwukicp+YTHSZK2OJADAMrl72fULzNBY+/prJevaiFJuuPjeer9InO2AQDlc6Jkz5LUwBiTaowJknSlpDEO5ACACjm5bFtL2QYAlM/VS/iNkjRdUkNjTJ4x5kZrbYmkOySNk7RM0qfW2iWuzAEAlaG8sv3N/M2ssw0AOM7Vq4tcdZrz30v63pWfGwBc5deyfWHT2vp+0Va9PGmV7h49X8+NX6lbu9TXpS0TFRzg73RMAICD3H7HRwBwV8dHtu/urOFDsxRVLVAPfrlIXZ7K1ds/r9OBIyVORwQAOMSrSrYxpp8xZnhhYaHTUQD4ED+/Y+tsf317B310YxvVjQnVv75dqo5PTtYrk1ap8GCx0xEBAFXMq0q2tTbHWjssMjLS6SgAfNCvO0iOHtZOX9zWTs2To/TMDyvV8YlJemrscu0oOux0RABAFXFinWwA8HpZKdF657poLdlSqNdy1+j1H9fonanrdGWrOhrWuZ4Soqo5HREA4EKUbABwoSYJkXp1cEut2V6k13PX6KMZGzRy5gZd2iJJw7rUU/1a7CYHAN7Iq6aLAIC7ql8rXM9cnqncP3XV4NZ19PX8zTr/uR918wezNWfDLqfjAQAqGSPZAFCFkmqE6tEBGbqzRwN9MG29PpixQeOX5isrpYZu6VxP5zeOc3zLdgDAuWMkGwAcEBMerPt6NtS0v3TXI/3Slb/3kIZ9OEfnP/+jRv+yUYeKjzodEQBwDijZAOCg0KAAXdchVbn3d9VLV7VQaJC//vLlInV8crJenbxahQdY/g8APJFXlWzWyQbgqQL8/dQ/M0E5d3TUyJvaKD2hup4et0Ltnpiof+Ys1eY9B52OCAA4A141J9tamyMpJzs7+2answDA2TDGqENajDqkxWjplr0aMWWt3p++Xu9PX69+zWrrpk71lJHIXgAA4O68qmQDgDdJT6iu5wc11/29Guqdn9dp9C8b9fX8LWqTGq02USXqVGrlz02SAOCWvGq6CAB4o8Soavp733RNe7CH/nZhY+XtPqiX5h1W92dz9e7UdSo6XOJ0RADASSjZAOAhIqsF6ubO9fTjn7rqD82DVTMsSI/mLFW7xyfq398tVd7uA05HBACUYboIAHiYAH8/tY4P0ANXdtC8jbv1ztT1emfqer398zr1yaitGzqmKiulhtMxAcCnUbIBwIO1qFNDL9epoQf7NNL709dr1MyN+m7RVjVPjtINHVPVJyNegf78oyUAVDVKNgB4gYSoanqwT2Pd1b2Bvpibp3enrtddo+apdmSIrm6boqta11F0WJDTMQHAZzC8AQBeJCw4QNe0q6uJ93XR29dmKzUmTE+PW6G2j0/U/Z8t0KI89hEAgKrgVSPZxph+kvqlpaU5HQUAHOXnZ9SjcZx6NI7Tqvx9en/6en05d7M+n5OnlnWidG37uuqTUVtBAYy1AIAreNV3V2ttjrV2WGQkGzUAwK8axEXosYubasZfe+gffdO1a/8R3T16vjo8OUnPj1+p/L2HnI4IAF7Hq0ayAQCnVz0kUDd0TNV17evqp1Xb9f609Xpp0iq9Onm1+jStrWvbpSgrpYaMYYMbADhXlGwA8DF+fkZdG8aqa8NYrd+xXx/O2KBPZ29SzoItapJQXde2r6v+mQkKCfR3OioAeCyvmi4CADgzdWPC9Pe+6ZrxYA/9+5IMFR8t1QOfL1Sb/5uof327VGu2FzkdEQA8EiPZAACFBQdoSJsUDW5dRzPW7tLImRv0wfRjG9y0q1dTV7dN0QXpcdwoCQAVRMkGABxnjFG7+jXVrn5Nbd93WJ/O3qSPZ27U7R/PVUx4sAa1StKVreooOTrU6agA4NYo2QCAU6oVEazbu6Xp1i719dOq7Ro5Y4Nez12j13LXqFvDWA1pU0ddG8bK348bJQHgZJRsAEC5/P2MujWMVbeGsdq856A++WWjRs/apBvfn63EqGq6qnWyrmiVrNiIEKejAoDbYHIdAKDCEqOq6b6eDTX1L931+pCWSo0J0zM/rFT7xyfpDyPnaOrqHbLWOh0TABznVSPZ7PgIAFUj0N9PfZrWVp+mtbVux359PHODPpuTp+8XbVPdmqEa1KqOBmYlMroNwGd51Ug2Oz4CQNVLjQnT3y46tgzg84MyFVs9RE+OXa52j0/SsA9ma9LyfJUcLXU6JgBUKa8ayQYAOCck0F+XtEjSJS2StGZ7kT6dvUlfzMnTD0vzFV89RJdnJ+mK7GRWJgHgEyjZAIBKV79WuB7s01j392yoicsK9MmsjXp18mq9PGm1OqbFaFCrZPVsEqfgAHaVBOCdKNkAAJcJ9PdT74x49c6I15Y9B/X5nDx9MmuT7hw1T1Ghgbq0RZIGtUpWw/gIp6MCQKWiZAMAqkRCVDXd1aOB7uiWpqlrdmj0rE36cMZ6vTN1nVrUidKVrZLVt1mCwoL50QTA8/GdDABQpfz8jDo1qKVODWppZ9FhfTVvs0bP2qQ/f7FIj+Ys1YVNa+uyrCS1rhstPza6AeChKNkAAMfUDA/WTZ3q6caOqZq7cY8+nbVJ3y3aqs/n5KlOdKgGtkzSpS0TuVkSgMehZAMAHGeMUVZKDWWl1NDD/dM1bsk2fT4nTy9MXKnnJ6xUu3o1dVlWkvo0jVdoED+6ALg/vlMBANxKaFDA8aUA83Yf0FdzN+vzuXn642cL9I9vFv93OklqtIxhOgkA90TJBgC4raQaobqzRwPd0T1Nszfs1uez8/Ttwi367ITpJAOzEpVUg+kkANwLJRsA4PaMMWpVN1qt6kb/z3SS5yccm07Svn5NDWyZpN4Z8axOAsAteNV3ImNMP0n90tLSnI4CAHCRk6eTfDl3sz6fc2w6yUNfL1bvjHhd3CJRHerXVIC/n9NxAfgoryrZ1tocSTnZ2dk3O50FAOB6STVCdVePBrqzbDrJl3M367uFW/TVvM2qFRGs/pkJuqRFopokVGf+NoAq5VUlGwDgm06cTvJI/3RNXl6gr+Zt1gfT1+vtn9epQWy4Lm6RqItbJCoxqprTcQH4AEo2AMCrBAf4q3dGbfXOqK09B47ou0Vb9fW8zXp63Ao9PW6F2qRG69KWieqdUVuR1QKdjgvAS1GyAQBeKyo0SEPapGhImxRt3HlA38zfrK/mbdafv1ikv3+zRBc0jtPFLRLV5bxaCgpg/jaAykPJBgD4hDo1/7sc4MK8Qn01b7NyFmzRd4u2qkZooC5sWlv9MxPUiu3cAVQCSjYAwKcYY5SZHKXM5Cj97aLGmrJqu76at0VfzM3TyJkbVTsyRH2b1Vb/zERlJHLDJICzQ8kGAPisQH8/dW8Up+6N4rT/cIkmLMvXmPlb9N609RoxZZ1SY8LULzNB/TMTlBYb7nRcAB6Ekg0AgKSw4AANaJ6oAc0TtefAEY1dvE1jFmzRy5NW6aWJq5Reu7r6N09Q32a12WESwO+iZAMAcJKo0CBd2bqOrmxdRwV7D+nbhVs1ZsEWPfGf5XriP8uVlVJD/TMTdGHT2qoVEex0XABuiJINAEA5YquH6IaOqbqhY6o27jygnIVblLNgix4es0SP5ixRh7QY9ctMUK8m8SwJCOA4SjYAABVUp2aobu+Wptu7pWll/j6Nmb9FYxZs0QOfL9RDXy1W14a1dFGz2urROE7hwfyIBXwZ3wEAADgL58VF6P5eDfXHnudpQV6hxszfou8WbdEPS/MVFOCnrudRuAFfxlc9AADnwBij5slRap4cpYcuaqw5G3fru4Vb9Z/FW/+ncNcLKFH24RIKN+Aj+EoHAKCS+PkZtaobrVZ1o/WPvun/W7j3HtY7S8czwg34CK/66jbG9JPULy0tzekoAAAfd3LhfuubSdriH/+bEW4KN+CdvOor2lqbIyknOzv7ZqezAADwKz8/o/Nq+GtY1ya/HeGmcANeia9iAACqULlTSpbmKzjAT53Pq6XeTeJ1fuM4RYayLCDgiSjZAAA45HSFe9ySbRq/NF8Bfkbt6tdUrybx6tkkTrERIU5HBlBBlGwAANzAiYX74X7pWphXqLFLtmns4m166OvF+vs3i5VVp4Z6Z8SrV5N4JUeztTvgzijZAAC4GWOMMpOjlJkcpQd6NdTK/CKNXbxNY5ds02PfLdNj3y1Tk4Tq6t0kXr0z4pUWGy5jjNOxAZyAkg0AgBszxqhhfIQaxkfo7vMbaMPO/RpXNsL97PiVenb8StWrFXa8cDdNjKRwA26Akg0AgAdJqRmmYZ3ra1jn+srfe0g/LDk2wv3mT2v1Wu4aJUZVU88mcerVJF7ZKTUU4O/ndGTAJ1GyAQDwUHHVQzS0XV0NbVdXu/cf0YRl+Rq3ZJtGztyod6euV43QQHVrFKue6XHq1KCWwlgaEKgyfLUBAOAFaoQF6fLsZF2enayiwyX6ccV2jV+6TROW5uvLuZsVFOCnjmkxuiA9Tj0ax7JSCeBilGwAALxMeHCALmpWWxc1q63io6WatW6Xxi/L1/il+Zq0vEDGSM2To3RBepwuaBzHjZOAC1CyAQDwYoH+fmqfFqP2aTH6R990Ld+2T+OXHivcT41doafGrlDdmqHHCnd6vLJSasjfj8INnCtKNgAAPsIYo8a1q6tx7eq6q0cDbS08qAlL8zV+WYHem7ZeI6asU3RYkLo3itUF6XHq1CBGoUFUBeBs8JUDAICPqh1Z7fiNk/sOFevHlds1fumxmyc/n5On4BPmcXdnHjdwRijZAABAESGB6tssQX2bJaj4aKl+Wbfr+LSSicsLJEmZSZHq3ihO3RvFqklCdfkxrQQ4LUo2AAD4H4H+fuqQFqMOaTF6uF+6lm3dp0nLj5XtFyau1PMTVio2IljdG8Wqe6NYdUiLYXlA4CR8RQAAgNMyxig9obrSE6rrju4NtLPosHJXbNek5QX6buFWjZ61SUH+fmpbv6Z6lJXu5OhQp2MDjqNkAwCACqsZHqyBWUkamJWkIyWlmr1+lyYtL9Ck5QV6eMwSPTxmiRrEhqt741j1aBSnlnWi2HUSPomSDQAAzkpQwH+XB3yob7rWbi/SpOUFmryiQG9PWac3f1yryP9v706DrKzuPI5//73QDXY3vSJIeqHpFlRUlgYEFElqTIxjtppsVZlkksykKslMTVIxM5XMi2yVmlQmyWSfSSYzZqksmsUsRo0LSlBEZXGDKDSrG4ZuNtHIfubFfbrtdBAavPS9jd9PVRf33uf0uedyOPDj9HnOGV3OoilNvGrqOC49u4naMaMK3WxpWJxWITsiXge8rqOjo9BNkSTpZae9qYr2pir+4ZJ2ntl3kLu6e1n8yHaWrNvOrx94ipKArtZ6Xjl1HIumNDF1fLWH4Oi0dVqF7JTS9cD1XV1d7yt0WyRJejmrqSznivMncMX5EzhyJPHgE7u5/dHtLH5kO5//3aN8/nePMr6mkkvPbmLRlCYWdDZSU1le6GZLeXNahWxJklR8SkqCGS11zGip46pXT+HpPfv4/frt/H59Dzeu2ca1Kx+nrCSY2VrXH7rPnVDjLLdGNEO2JEkaVuPHVvK22S28bXYLBw8f4f7HdrNk3XaWrOvhCzev4ws3r2NcdUUWuMdxcUcjY8c4y62RxZAtSZIKpry0hDmT6pkzqZ5/vXwq25/Zx+/X97BkfQ83r32an616gtKSYEZzLYum5EL3uRM8CEfFz5AtSZKKxriaSt7S1cxbupo5dPgIDzy+Oxe61/XwxVvW88Vb1tNYVcHCsxtZNGUcCzsb3bFERcmQLUmSilJZaQldbfV0tdVz1aun0LN3P0uzWe7bH93OdaufpCTgwuZaFnY2sfDsRi58hftyqzgYsiVJ0ojQVP3CQTiHsx1LlqzrYen6Hr5+ezdfXdxNdUUZ8zsauKSziYWdTbQ0ePqkCsOQLUmSRpzSkmBmSx0zW+r4yGVns/tPB7h74w7u7O5h6fpebl77RwBaG8ZwSWcjl3Q2MW9yg9sEatgYsiVJ0ohXO2ZU/77cKSU29z7Hnd293Nndwy9XP8kP73ms/wbKSzqbuOTsRi6YONalJTplDNmSJOm0EhH9p0/+3fw2Dhw6wv2P7eoP3V9ZvJ4v37aemsoyFnTkZrkv6Wykud6lJcofQ7YkSTqtjSorYW57A3PbG/joa6aw67kDLNvYy53re1na3cNNa54GoL3xjP6lJRdNbihwqzXSGbIlSdLLSt0Zo7jygrO48oKzSCmxsec57uzu4c7uXn668gm+v3wrZSVB+9jggUPrWdDRyPTmWspdWqITYMiWJEkvWxFBx7gqOsZV8Z4Fk9h/6DCrt+7mzu4ebrp/M19b3M1XbutmzKhS5kyq5+KORuZPbmTq+GoPxNExGbIlSZIyFWWlzJvcwLzJDcypfJoZcxawfNMO7t7Yy7INvXz2hkcAqD9jFPMmN3BxRyMLJje6VaD+giFbkiTpRYwdU87l08Zz+bTxADy9Zx/LNvSyLAvdNzy0DYBXHrklMwAADaxJREFU1I1mweRGFnQ2Mn9yA41VFYVstoqAIVuSJGmIxo+t7D8Qp289990be7mru5cb12zj2pWPAzB1fDXzJzdycWcDcyY1UFVh5Hq5scclSZJOwsD13O+a18bhI4k1T+7pn+X+4b1buXrZZspKgguba1kwuSF3E2VLLRVlpYVuvk4xQ7YkSVIelGZh+sLmWj64qIN9Bw+zeusu7trQy7KNO/jGHRv42u0bqCwvYVZrHfPac2u/z59Yy6gydy453RiyJUmSToHK8lLmdzQyv6MRgD3PH+SeTTtYvnEH92zawRdvWQ/A6PJSutrqcjdctjdwvidRnhYM2ZIkScNg7OhyXnPeeF5zXu4myp3PHeDeTbnAvXzTDv7jd+sAOGNUKbMn1TOvvYGL2hs476waQ/cIZMiWJEkqgPozRvHa8yfw2vMnAND77H7u6QvdG3ewZF0PANUVZcyZVM9F2fKScybUUOoe3UXPkC1JklQEGqsq+k+iBNj+zD7u2byzf3nJ4ke3A1BTWcbcbJZ7XnuDB+MUKUO2JElSERpXU8nrLzyL11+YC91P79n3wpruzTu49Q9/BKB2TDlzJ9Uzd1IDc9vrmTreme5iYMiWJEkaAcaPreSNMybyxhkTAXhy9/Pcs/GFNd03r82F7urKMma31TNnUj1zJ9UzbeJYyl3TPewM2ZIkSSPQxNrR/QfjQC50r9i8k3s37+DezTu5PVteMrq8lFmtdf2h+8LmWirL3af7VDNkS5IknQYm1o5m4oCZ7p69+1mxZSf3bsqF7i/ftp6UYFRpCdOba5nbnpvtntlSxxmeSJl3p9XvaES8DnhdR0dHoZsiSZJUUE3VFVxx/gSuyHYv2f2nA6zYsov7Nu/gvs07+a8lG/n67RsoLQmmTRzLRZNyoburrZ6xo8sL3PqR77QK2Sml64Hru7q63lfotkiSJBWT2jGjuOzcM7ns3DMBeHb/IVZtfSF0f3fZFr69dBMRMHV8TXYzZT2zJ9XTWFVR4NaPPKdVyJYkSdLQVFWUcenZTVx6dhMA+w4e5v7HdnPf5p3ct2UH16x4jO/dvQWAyU1n5Ga5W+uZ3VZPc/1oItzB5FgM2ZIkSaKyvDR3tPvkBqCTA4eO8PCTe3Khe/MOfvvQNn5y3+MAjKuuYHZbPbNa65jdVs85E6o9lXIQQ7YkSZL+wqiyEma11jGrtY4PLJrMkSOJ9dv3smLLLlZu2cnKLbu44eFtAIwZVcrMljq62nKhe3pz7cv+ZsqX96eXJEnSkJSUBFPH1zB1fA3vvKgVgKd2P8/KrbnQvWLLLr66uJuUoLQkOO+smmx5SR2z2uoYV11Z4E8wvAzZkiRJOiln1Y7m9bWj+0+lfGbfQVZv3cWqrbtYsWUnP75vK1cv2wxAa8OY/tDd1VbP5KYzTut13YZsSZIk5UVNZTmLpoxj0ZRxABw4dIS1T+1h5ZZc6L5j3XZ+sfoJAOrGlNPVls10t9YzbWINFWWnzyE5hmxJkiSdEqPKSpjRUseMljret7CdlBKbe5/rD90rt+7i1j/8sb/sBRPHMqu1jpmtdcxsqaOpeuRuHWjIliRJ0rCICNqbqmhvquKts5uB3MmUq7buZFW2zKRvv27ILTGZ1ZJb0z2rtY7OcdWUloyMJSaGbEmSJBVMU3UFl0+bwOXTcidT7jt4mLVP7ekP3Uu7e7nu/icBqK4oY3pLbf+uJ9Oba6muLM7TKQ3ZkiRJKhqV5aXMaq1nVms9ACklHt/5PKse65vt3s3XFndzJEEETDmzmrd2NfPeiycVuOV/zpAtSZKkohURtDSMoaVhDG+a8QoA9u47yIOPZ7Pdj+3i+YOHC9zKv2TIliRJ0ohSXVnOxZ2NXNzZWOimvCjPv5QkSZLyzJAtSZIk5ZkhW5IkScozQ7YkSZKUZ4ZsSZIkKc8M2ZIkSVKeGbIlSZKkPDNkS5IkSXlmyJYkSZLyzJAtSZIk5ZkhW5IkScozQ7YkSZKUZ4ZsSZIkKc8M2ZIkSVKeGbIlSZKkPDNkS5IkSXlmyJYkSZLyzJAtSZIk5VmklArdhryLiB5ga6Hb8TLWCPQWuhE6IfbZyGOfjTz22chjn408w91nrSmlpqNdOC1DtgorIlamlLoK3Q4NnX028thnI499NvLYZyNPMfWZy0UkSZKkPDNkS5IkSXlmyNap8D+FboBOmH028thnI499NvLYZyNP0fSZa7IlSZKkPHMmW5IkScozQ7ZOSkRcHhHrImJDRHzsKNcXRcSeiHgg+/pEIdqpF0TE1RGxPSLWvMj1iIivZX36UETMHO426s8Noc8cZ0UmIpoj4o6IeCQi1kbEh45SxrFWRIbYZ461IhIRlRFxX0Q8mPXZp49SpuDjrGy431AjX0SUAt8ELgOeAFZExG9SSn8YVPTOlNKVw95AvZjvAd8AfvAi118LdGZfc4H/zn5V4XyPY/cZOM6KzSHgqpTS6oioBlZFxK2D/n50rBWXofQZONaKyX7gVSmlZyOiHLgrIm5KKd0zoEzBx5kz2ToZc4ANKaVNKaUDwDXAGwrcJh1HSmkpsPMYRd4A/CDl3APURsSE4WmdjmYIfaYik1LallJanT3eCzwCTBxUzLFWRIbYZyoi2dh5Nntann0Nvsmw4OPMkK2TMRF4fMDzJzj6X0jzsh/l3BQR5w1P0/QSDLVfVVwcZ0UqItqAGcC9gy451orUMfoMHGtFJSJKI+IBYDtwa0qp6MaZy0V0MuIorw3+H+RqckeNPhsRVwC/IvcjGxWvofSriovjrEhFRBXwC+DDKaVnBl8+yrc41grsOH3mWCsyKaXDwPSIqAV+GRHTUkoD718p+DhzJlsn4wmgecDzVwBPDSyQUnqm70c5KaUbgfKIaBy+JuokHLdfVVwcZ8UpWyP6C+BHKaXrjlLEsVZkjtdnjrXilVLaDSwBLh90qeDjzJCtk7EC6IyISRExCng78JuBBSJifERE9ngOuT9rO4a9pToRvwHeld2RfRGwJ6W0rdCN0otznBWfrD/+D3gkpfSfL1LMsVZEhtJnjrXiEhFN2Qw2ETEa+Cvg0UHFCj7OXC6iE5ZSOhQR/wTcDJQCV6eU1kbE+7Pr3wLeDHwgIg4BzwNvT558VFAR8RNgEdAYEU8AnyR3s0hfn90IXAFsAP4EvKcwLVWfIfSZ46z4LADeCTycrRcF+DegBRxrRWoofeZYKy4TgO9nu52VAD9NKf12UA4p+DjzxEdJkiQpz1wuIkmSJOWZIVuSJEnKM0O2JEmSlGeGbEmSJCnPDNmSJElSnhmyJekUiYgvR8SHBzy/OSL+d8DzL0XER/L8ns/ms76szunZKXd9zz8VER8dwvdFRNweETV5aMOoiFgaEW49K2lEMGRL0qlzNzAfICJKgEbgvAHX5wPLCtCuEzWd3H6zJ+oK4MGjHFF9wlJKB4DFwNteal2SNBwM2ZJ06iwjC9nkwvUaYG9E1EVEBXAOcH9EVEXE4ohYHREPR8QbACLi8xHxwb7Kshnkq7LH/xIRKyLioYj49NHe/GhlIqItIh6JiO9ExNqIuCU7MY2ImJ2VXR4RX4iINdmprp8B3hYRD0REX8g9NyKWRMSmiPjnF/n87wB+PeB91wxo20cj4lPZ4yXZrP/SrG2zI+K6iOiOiM8OqO9XWZ2SVPQM2ZJ0iqSUngIORUQLubC9HLgXmAd0AQ9lM7T7gDellGYCrwS+lB3hfA1/PnP7VuBnEfFqoBOYQ26WeVZELBz43scp0wl8M6V0HrAb+Jvs9e8C708pzQMOZ5/hAPAJ4NqU0vSU0rVZ2anAa7L6PxkR5Uf5LVgArBrib9eBlNJC4Fvkgvk/AtOAd0dEQ1ZmDTB7iPVJUkEZsiXp1Oqbze4L2csHPL87KxPAv0fEQ8BtwETgzJTS/cC4iDgrIi4EdqWUHgNenX3dD6wmF3g7B73vscpsTin1HR+9CmiLiFqgOqXU16YfH+dz3ZBS2p9S6gW2A2cepUx9Smnvcerp85vs14eBtSmlbSml/cAmoBkgpXQYOBAR1UOsU5IKxhtIJOnU6luXfT65mdjHgauAZ4CrszLvAJqAWSmlgxGxBajMrv0ceDMwntzMNuRC+edSSt8+xvsetUxEtAH7B7x0GBidlT8Rg+s42r8nhyKiJKV0BDjEn0/sVA4q21ffkUF1HxlUdwW5mX9JKmrOZEvSqbUMuBLYmVI6nFLaCdSSWzKyPCszFtieBexXAq0Dvv8a4O3kgvbPs9duBt4bEVUAETExIsYNet+hlOmXUtpFbr34RdlLbx9weS9wMrPH64D27PEfyc3KN2Tr0a880cqyZSM9KaWDJ9EWSRpWzmRL0qn1MLldRX486LWqbKkFwI+A6yNiJfAA8GhfwZTS2mx5xJMppW3Za7dExDnA8tzSbZ4F/pbcsg2OU+bwMdr698B3IuI5YAmwJ3v9DuBjEfEA8LkT+Ow3AIuADdl/ID5Dbk365oGf8QS8ErjxJL5PkoZdpJQK3QZJUhGIiKqU0rPZ448BE1JKH3oJ9U0AfpBSuixP7bsO+HhKaV0+6pOkU8mZbElSn7+OiI+T+7dhK/Dul1JZSmlbtlVgzUvdKzvbSvBXBmxJI4Uz2ZIkSVKeeeOjJEmSlGeGbEmSJCnPDNmSJElSnhmyJUmSpDwzZEuSJEl5ZsiWJEmS8uz/AZIKfOP5Uy/dAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure(figsize=(12,8))\n",
    "semilogy(wavelength, specAOD)\n",
    "grid()\n",
    "xlabel('Wavelength (um)')\n",
    "ylabel('Spectral Optical Depth')\n",
    "legend(['Layer 1', 'Layer 2', 'Layer 3', 'Layer 4'], loc='best')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The spectrally-integrated aerosol optical depth is   1.00 \n"
     ]
    }
   ],
   "source": [
    "AOD = trapz(specAOD, wavelength)\n",
    "print('The spectrally-integrated aerosol optical depth is %6.2f ' % AOD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "file_extension": ".py",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  },
  "mimetype": "text/x-python",
  "name": "python",
  "npconvert_exporter": "python",
  "pygments_lexer": "ipython3",
  "version": 3
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
